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We combine density functional theory within the local density approximation, the 
quasiharmonic approximation, and vibrational density of states to calculate single 
crystal elastic constants, and bulk and shear moduli of diamond at simultaneous 
high pressures and temperatures in the ranges of — 500 GPa and — 4800 K. Com- 
parison with experimental values at ambient pressure and high temperature shows 
an excellent agreement for the first time with our first-principles results validating 
our method. We show that the anisotropy factor of diamond increases to 40% at high 
pressures and becomes temperature independent. 
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Diamond has singular physical properties. Its strength and transparency make it re- 
markably useful for industrial applications in abrasive machining, micro electromechanical 
devices etc. 1 It is also an anvil material for experimental studies under extreme pressures 
and temperatures. Pressures greater to that reached at the center of the Earth (~ 364 
GPa) can be produced in diamond anvil cells (DACs) using single crystal (SC) diamonds as 
opposing anvils. Applications of diamond to geoscience and planetary science are ideal to 
study high pressure-temperature behavior of magmas and minerals relevant to the Earth and 
other planets. For example, in 2004 a high pressure form of (Fe,Mg)Si03, postperovskite, 
was discovered using DAC techniques 2 , this mineral exists above the mantle-core boundary 
of the Earth. 

Natural diamond can be SC or polycrystalline (PC). The high quality SC diamond is 
optically transparent and often used as gems. It also possesses higher strength than the 
PC type, which is more opaque or translucent with impurities and inclusions. PC diamond, 
however, is ideal for cutting tools due to its toughness. 

Thin films of PC diamond with very low growth rate (~ 1/im/h) were synthesized and 
reported in the 1960s using chemical vapor deposition (CVD), but those films were of poor 
quality. SC diamond synthesis is performed near graphite-diamond phase boundary condi- 
tions (~ 5 — 6 GPa and > 1500 K) within large- volume presses 3 . The main limitations to 
this production are the size of the press, yielding only diamonds of some millimeters and 
a very slow growth rate 4 . However, in recent years SCs of diamonds larger than 2 cm and 
more than 10 carats have been successfully fabricated using microwave plasma CVD tech- 
niques. The purest CVD diamond was reported to have comparable features to those of high 
purity natural SC with hardness and toughness that can be tuned by controlling growth and 
annealing 5 ' 6 . In 2003 in Japan, using large- volume presses 7 , graphite was directly converted 
into diamond to form nanopolycrystalline diamond (NPD) at simultaneous high pressure 
and temperature (~ 15 GPa and > 2600 K) within minutes. NPD was found to have high 
elastic stiffness, very high fracture toughness and similar or higher hardness than that of 
most single crystal diamonds. These features and the fact that NPD rods can be fabricated 
up to 1 cm in diameter and length (14.5 carats), make NPD very suitable for larger anvils 
than the usually available in DACs. Larger sample volume in high pressure apparatuses 
has open new windows into mineral physics studies of materials that require accurate and 
simultaneous very high pressures and temperatures such as those found deep in the Earth 
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and other planets 8 . 

With these new techniques to fabricate diamonds and their increasing applications, ac- 
curate elasticity of diamond at simultaneous high pressures and temperatures acquires new 
significance. However, it cannot be found in the literature in neither experimental form 
nor computational modeling. Using Brillouin scattering, elastic moduli of diamond were 
measured as a function of temperature in the range of 300-1600 K at ambient pressure 9 . 
A molecular dynamics study attempted the computation of elastic constants between 100 
and 1100 K with a modest degree of success 10 . Single crystal elasticity and bulk modulus 
of diamond were computed as a function of pressure at 300 K in the range of 0-500 GPa 11 , 
using first-principles methods based on density functional theory (DFT) 12 ' 13 as implemented 
in the CRYSTAL03 (CR03) program 14 and a qusiharmonic Debye model (QHDM) 15 . But 
the lack of experimental and other computational values for such quantities made impossible 
a comparison for validation of these CRYSTAL03 results. 

In this letter, for the first time we address the elasticity of single crystal diamond and 
of a polycrystalline aggregate at simultaneous high pressures and temperatures. The com- 
putational approach is based on DFT 12,13 within the local density approximation (LDA) 16 . 
At arbitrary pressures, cubic equilibrium structures of diamond (space group Fd3m, two C 
atoms/primitive cell) were found using the variable cell-shape damped molecular dynamics 
approach 17 ' 18 as implemented in the quantum-ESPRESSO (QE) code 19 . The C pseudopo- 
tential was of Vanderbilt-type with cutoff radii of 1.3 a.u., and partial core correction. We 
used a plane-wave kinetic energy cutoff of 50 Ry and 200 Ry for the charge density. The 
k-point sampling for charge density was determined on a 4 x 4 x 4 Monkhorst-Pack grid 
of the Brillouin Zone (BZ) shifted by (1/2,1/2,1/2) from the origin. These parameters 
corresponded to having interatomic forces and stresses smaller than 10~ 4 Ry/a.u. 

To obtain the independent cubic static elastic constants Cu, C 12 , and C 44 (in Voigt 
notation 20 ) at each pressure, the equilibrium structures were strained, and their internal 
degrees of freedom re-relaxed. Then elastic constants were extracted using the relationship 
Cjj = CijkiSki- Positive and negative strains of 1% magnitude were applied in order to attain 
accuracy in the limit of zero strain. This method has been applied successfully for over a 
decade 21-24 . 

Dynamical matrices were calculated using density functional perturbation theory (DFPT) 5 
At each pressure, a dynamical matrix was obtained on a 4 x 4 x 4 q-point mesh. Force 
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FIG. 1. (Color online) Temperature dependence of (a) Cn; (b) C12; (c) C44; and (d) bulk, K, and 
shear, G, moduli. Dashed lines indicate a lesser degree of confidence within the QHA limit. 



constants were extracted and interpolated on a 12 x 12 x 12 mesh to produce VDoS. 

Then we exploited the dependence of phonon frequencies on compression to determine 
the thermal contribution to the Helmholtz free energy F, within the QHA 20 , 



F (e, V, T) = U st (e, V) + ^ q , m (e, V) + 



q,m 



q,m 



exp 



k B T 



(1) 



TABLE I. Comparison of elastic mouduli of diamond at ambient conditions (P=0 GPa and T= 
300K) 



Reference C n (GPa) C 12 (GPa) C44 (GPa) K VRH (GPa) G VRH (GPa) 

This work LDA 1078.99 140.81 577.37 453.54 531.31 

GGA 26 1060 125 562 436.67 522.08 

MD 10 963.6-982.2 148.4-165.8 674.5-694.0 425.76-437.93 550.12-565.99 

B3PW 11 1097.5 115.5 598.2 442.8 552.74 

Exp 27 1079 124 578 442 535 

Exp 28 1076.4 ± 0.2 125.2 ± 2.3 577.4 ± 1.4 442.3 534.27 

Exp 9 1080.4 127.0 576.6 444.8 534.32 



where q is the phonon wave vector, m is the normal mode index, T is temperature, U s t 
is the static internal energy at equilibrium volume V under isotropic pressure P and in- 
finitesimal strain e, h and k R are Planck and Boltzmann constants, respectively. Isother- 
mal elastic constants are given by Cj^ u = [d 2 G(P,T) /deijdeki] PT with G = F + PV, the 
Gibbs energy and i, j, k, I = 1, . . . , 3. To convert to adiabatic elastic constants, we use 
Cijki = Cjj kl + (T/VCv) (dS/deij) (dS/deu) SijSu, (Cy =heat capacity at constant V, and 
S =entropy). Cfj kl have the static-CV,^ and phonon-CV^ contributions, the latter can be 
expressed as function of strain and mode Griinesein parameters. Assuming that the angular 
distribution of phonons is isotropic, isothermal elastic constants can be calculated without 
performing phonon calculations for strained configurations 29 . This approximation is equiva- 
lent to assuming that thermal pressure is isotropic. Even though this assumption is not com- 
pletely rigorous 30 , the method is quite accurate compared to experimental uncertainties 31 . 
Bulk and shear moduli were calculated as Voigt-Reuss-Hill averages 32 using the adiabatic 
Cy, i.e., K VRH = (K v +K R )/2 and G VRH = (G v +G R )/2, whereK v = K R = (C n + 2C 12 )/3, 
G v = (C n - C 12 + C 44 )/5, and G R = 15 [12/(C n - C 12 ) + 9/C u )}-\ 

The computed elastic constants Cn, C12, and C44, and bulk Ky R H and shear Gvrh, 
moduli as a function of temperature are shown in Figure 1. We compared our results with 
experimental trends in the temperature range of 300 — 1700 K at ambient pressure and 
we found excellent agreement. The previous molecular dynamics calculation 10 was quite 
unsuccessful at reproducing the experimental values for Cij individually. The quality of our 
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findings with respect to measurements validates our overall computational approach. 

Then we obtained the elasticity as a function of pressure between and 500 GPa at 
several temperatures, (Figure 2). Unfortunately experimental values are only available at 
ambient conditions for comparison. The only other calculation of elasticity as a function 
of pressure 11 reported using CR03 and the QHDM was compared only to measurements at 
room temperature. Figure 2(a) shows that CR03 results are quite different from ours. There 
is agreement in the bulk modulus trend [Figure 2(b)], but that probably reflects the fact 
that the compression curves are similar in both studies. 

Since the QHA is valid up to about two thirds of the melting temperature (for diamond 
Tmeit = 3823 K at P = GPa and it increases in the range — 500 GPa 33 ), dashed-lines in 
Figures 1 and 2 indicate a lesser degree of confidence for the QHA results 34 . Table I compares 
results on the elasticity of diamond at ambient conditions with the other DFT-based values 
and experimental measurements. One can see clearly the improvement in the computation 
of Cij, Kvrh, and Gvrh with respect to the other computational methods. Our computed 
bulk modulus is 2.5% larger than experimental values. This is typical of LDA results that 
underestimate volume by ~ 1% and reflects on C\i only. Finally, we solved the Christoffel 
equation 35 for acoustic velocities in single crystal, det \Cijkifijni — pV 2 5ik\ = 0, where V is 
elastic wave velocity, n is propagation direction, 8^ is Kronecker delta, and p is density. 
For a given n, there are three solutions, i.e., one P-wave (Vp) and two S- waves (Vsi, V52). 
Figures 3(a-c) show the variation with crystallographic direction (anisotropy) of Vp, Vsi and 
V52 at and 100 GPa and 300 and 1000 K for comparison. We found similar qualitative 
behavior independent of pressure and temperature for the three waves. In Figure 3(d) we 
show the anisotropic factor of diamond A [(2C44 + C\<i)jC\\ — 1] x 100, as a function of P 
at several Ts. One can observe that A becomes roughly independent of T and P. A often 
decreases with P 36 , but this is not the case for diamond. As seen in Figures 3(a-c), velocities 
vary with direction much more at high pressures. 

In conclusion, we have computed the high pressure and temperature elasticity of diamond 
by first principles using a novel analytical approach. In this new powerful method, vibra- 
tional density of states for strained configurations are unnecessary 29 . Treatment of strain 
Griineisen parameters via isotropic averages reduced greatly the computational cost of the 
task, which, otherwise, would have been much more demanding. We have shown an out- 
standing improvement in the calculated elasticity of diamond as a function of temperature 
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FIG. 2. (Color online) Pressure dependence of (a) Cn, C12, and C44; and (b) bulk, Kvrh, and 
shear, Gvrh, moduli. Dashed lines indicate a lesser degree of confidence within the QHA limit. 



at ambient pressure (Figure 1). These results increase the confidence in our approach to 
calculate high temperature elasticity (Figure 2). In general, accurate elasticity of diamond 
at high pressures and temperatures is a step forward in improving and opening applications 
of diamonds in science and technology, particularly in geophysics and planetary sciences 
where extreme conditions of pressure and temperature are typical. The elastic anisotropy 
of diamond is shown to increase considerably at high pressures. This result is significant 
for understanding the performance under pressure of large SC and PC diamonds specimens 
nowadays synthesized by CVD or as NPDs. 
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